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This paper presents results of a theoretical investigation of transport in a numerical model of a 
two-dimensional Kolmogorov flow. We investigate the changes in its mixing properties associated 
with transition from laminar regime to turbulence. It is found that significant changes in the flow 
do not always lead to comparable changes in its transport properties. On the other hand, some very 
subtle changes in the flow can dramatically alter the degree of mixing. We show that interaction of 
multiple resonances can provide an explanation for many of these seemingly paradoxical results. 



I. INTRODUCTION 

Two-dimensional (2D) flows have proven to be very 
useful for studying various phenomena in fluid dynam- 
ics since, compared to three-dimensional flows, they are 
much more amenable to theoretical analysis and numer- 
ical investigation. In particular, much of our fundamen- 
tal understanding of transport properties of fluid flows 
has been developed using 2D models. Effectively 2D 
fluid flows are responsible for transport and mixing in 
many geophysical processes such as atmospheric [TJ |2] 
and oceanic [3 , 4 flows as well as in convection processes 
within the Earth's mantle [5l|6]. Two-dimensional lami- 
nar mixing is a key process in many types of microfluidic 
essays [71 18] , such as ones used for gene expression profil- 
ing [9], and in numerous technological applications, such 
as the production of polymer blends [10 . The reduc- 
tion to two dimensions has also provided insights into 
many difficult 3D problems ranging from mixing in the 
radiation zones of rotating stars [11 to confinement of 
thermonuclear plasmas [12 . 

Much of our understanding of transport properties of 
fluid flows comes from experimental observations or nu- 
merical simulations of the advection, or stirring, of pas- 
sive tracers by the flow. The dynamics of passive tracers 
in 2D flows of incompressible fluids is formally described 
by a Hamiltonian system with one degree of freedom 

y = vy = -5^^, (1) 

where the stream function ^{x^y^t) plays the role of a 
Hamiltonian and the coordinates x and y are the conju- 
gate variables. 

Time-independent 2D flows are integrable (because ^ 
is conserved), and the trajectories of tracers coincide with 
the closed stream lines of the flow, resulting in poor 
mixing. The introduction of time-dependence effectively 
makes the velocity field a three-dimensional flow, which 
is generally nonintegrable. Stream lines in such flows can 
be chaotic even if the underlying velocity field is regular 
(e.g., stable and time-periodic). If this is the case, the 
stream lines will diverge exponentially fast from one an- 
other, resulting in rapid stretching and folding of fluid 
elements. This process, known as chaotic advection or 



Lagrangian chaos, is the underlying mechanism respon- 
sible for dramatically improved mixing. 

One of the simplest 2D models in which chaotic mixing 
can occur is the 'blinking-vortex' flow studied by Aref 
[13] . Historically the first study of the mixing properties 
of a fluid system, this model was originally proposed as 
an idealization of a periodically stirred fluid and consists 
of a pair of spatially separated fixed point vortices which 
are alternately turned on for one half of the period T. 
Numerical simulations showed that when both vortices 
are running continuously (i.e., T = 0), the flow is time- 
independent and, thus, integrable. For small nonzero T, 
it was found that the trajectories nearest the vortices 
become chaotic. The size of the mixed region increases 
monotonically with T until, at some finite critical value 
of T, the entire domain becomes uniformly mixed. 

The first analytic investigation of mixing in a time- 
periodic 2D flow is due to Khakhar et al [14 . This study 
introduced an idealized model, known as the 'tendril- 
whorl' flow, in which uniform shear is followed by differ- 
ential rotation and showed that mixing takes place in the 
vicinity of separatrices associated with saddle fixed points 
of the period-T map of the flow, while the KAM tori sur- 
rounding the elliptic fixed points serve as transport barri- 
ers. The same structures were shown to also control mix- 
ing in the 'blinking vortex' flow. These studies demon- 
strated that laminar, time-periodic, area-preserving 2D 
flows can produce efficient mixing. 

The results of these idealized models raised questions 
as to whether or not real-world laminar fluid flows could 
give rise to chaotic stream lines. This prompted the 
analytical and numerical study of chaotic advection in 
a journal bearing Stokes flow with physical boundary 
conditions [151 HI]- The basic setup is that of a Cou- 
ette flow between non-coaxial rotating cylinders, where 
time-periodicity is introduced by alternating the rota- 
tion between the inner and outer cylinders. By varying 
the distance between the axes of the cylinders as well as 
the time-interval for which one of the cylinders rotates, 
one can obtain various flow patterns with both regular 
and chaotic trajectories. The experimental realization of 
this flow [IT showed excellent agreement with numeri- 
cal results. Subsequently, the experimental study of cav- 
ity flows by Chien et al. [18^ showed the existence of 
transverse intersections of homo/heteroclinic manifolds 
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at small Reynolds numbers, providing more evidence for 
the mixing capabilities of 2D laminar flows. 

Rom-Kedar et al. [19 proved the existence of chaotic 
trajectories analytically for a model flow produced by a 
pair of time-independent point vortices perturbed by a 
time-periodic shear. The theory of lobe dynamics devel- 
oped in this paper set up a framework for quantitative 
description of transport across separatrices of the unper- 
turbed flow which evolve into a homoclinic tangle in the 
presence of perturbation. In conjunction with the ana- 
lytic techniques introduced by Melnikov [20 , this frame- 
work enabled them to estimate the fluxes between differ- 
ent regions of the flow domain. 

Recently, most experimental and many theoretical 
studies of mixing in 2D flows have used thin layers of 
electrolyte placed over various arrangements of perma- 
nent magnets. The fluid flow is driven by the Lorenz 
force which arises when electric current flows through 
the electrolyte. Since this setup is closely related to our 
work, we describe here other studies which used it. 

Rothstein et al. [21 discovered the existence of per- 
sistent spatial patterns, which they called strange eigen- 
modes, in a flow driven by a combination of time-periodic 
current and either a disordered or a square array of mag- 
nets. These patterns were shown to emerge as a result 
of a delicate balance between advective stretching and 
molecular diffusion. The process of mixing was observed 
to continue even after these structures reached an asymp- 
totic shape. The same experimental setup was subse- 
quently used to investigate the rate of mixing [22]. By 
examining the spatial structure of persistent patterns, 
it was found that locally, mixing rates are controlled by 
stretching, but on large scales they are governed by diffu- 
sive transport. Additionally, it was discovered that mix- 
ing rates could be dramatically enhanced by breaking 
certain spatial and temporal symmetries. 

Voth et al. [23] used a disordered array of magnets 
and time-periodic current to drive the flow and were 
able to use flow measurements to construct forward and 
backwards finite-time Lyapunov exponent (FTLE) fields 
which follow the time-evolution of the unstable and sta- 
ble manifolds of saddle points of the flow, thus provid- 
ing an empirical method for visualizing the geometrical 
structures underlying the mixing process. A follow-up 
experimental study carried out using magnets arranged 
in a square, hexagonal, and a disordered array [24 found 
that the probability distribution of FTLEs exhibited self- 
similar behavior regardless of the flow pattern or the de- 
gree of mixing in the system. 

Fluid mixing was also studied in time-dependent flows 
driven by steady current. Danilov et al. [25] performed a 
combined experimental and theoretical study of mixing 
by a time-periodic four- vortex flow. Numerical study of a 
truncated analytic model showed that separatrices parti- 
tioned the flow domain into regions with different mixing 
rates and that transport between these regions was a rel- 
atively slow process compared to the mixing within these 
regions. The theory of adiabatic chaos was used to ex- 



plain the results and show that long-term transport could 
effectively be modeled as a random walk of an adiabatic 
invariant. 

This paper investigates mixing properties of a range of 
2D flows arising as intermediate stages in the transition 
from the so-called Kolmogorov flow [26l [27] to turbu- 
lence. Unlike the majority of other studies of mixing in 
2D flows where the time-dependence is due to external 
monochromatic forcing, our focus is on time-dependence 
that arises naturally as a result of fluid-dynamic insta- 
bilities, producing flows ranging in their temporal com- 
plexity from time-periodic to quasi-periodic and chaotic. 
The paper is organized as follows: In Sect. [ll] we in- 
troduce the model of the fluid flow and characterize the 
flow states that emerge in the transition from the laminar 
to the turbulent regime. The mixing properties of these 
flows are described and analyzed in Sect. |III[ Finally, 
summary and conclusions are presented in Sect. |IV| 



II. PROBLEM DESCRIPTION 

A. Model of the Fluid Flow 

We consider a model of an experimental flow described 
in Ref. [28] which employs bar magnets with alternating 
polarity to generate a Kolmogorov flow in a layer of elec- 
trolyte supported by a liquid dielectric. The flow in the 
conducting layer can be described by the following equa- 
tion for the vorticity Q = — V^^: 

dtn + /3v • V^] = -an^A sin ky (2) 

where k = tt/w and w is the width of individual magnets. 
Parameters P = 1^ u = 0.0115 cm^/s, and a = 0.1141 
s~^ were selected to be representative of a typical ex- 
perimental setup. Furthermore, we chose the domain 
width Ly = 5 cm corresponding to four magnets of width 
w = 1.25 cm and the length = 2Ly = 10 cm. For sim- 
plicity, unlike the experimental system which is larger 
and features physical (no-slip) lateral boundary condi- 
tions, we assume periodic boundary conditions. The ef- 
fect of the bottom boundary, however, is included in our 
model via the Rayleigh friction term —aQ. The impor- 
tance of this term is described by the non-dimensional 
combination F = a/vk'^ ~ 1.57 which shows that it is 
comparable to the viscous term vV^^. Finally, A mea- 
sures the strength of the driving force and is used as a 
control parameter analogous to the Reynolds number. 

The vorticity equation ^ was solved numerically using 
a spectral (Fourier) method with 64 x 128 modes. As 
a check, we recalculated the bifurcation sequence using 
128 X 256 modes which yielded less than a 1% difference in 
both the leading stability eigenvalues and the location of 
the bifurcations. Temporal discretization used a second- 
order, implicit-explicit, operator-splitting scheme with an 
adaptive time step [29] . The Crank-Nicolson method was 
used for the linear and forcing terms in Fourier space. 
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FIG. 1. Laminar flow L at A = 0.1 (a) and spatially mod- 
ulated flow M at A = 0.250 s"^ (b). Velocity field (arrows) 
is overlayed on top of vorticity field (grayscale). 
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FIG. 2. Bifurcation diagram. The relative vorticity magni- 
tude Qo = ||Q — II2 — is shown, where c is a constant 
chosen to separate the various branches of the diagram for 
visualization purposes. Solid and dashed lines denote stable 
and unstable states, respectively. Periodic orbits are repre- 
sented by their time- averaged values. Inset shows the region 
where the P3 branch exists. 



while the fourth-order Runge-Kutta method was used for 
the advection term in real space. The use of the so-called 
Strang-Marchuk splitting [30 ensures that the resulting 
scheme is second order in time. 



B. Prom Kolmogorov Flow to Turbulence 

In this section we describe the transition to 2D tur- 
bulence in our model system as the value of the control 
parameter A is increased. Unlike many shear flows in 
3D which transition directly from laminar flow to turbu- 
lence, here we find a rather complicated sequence of tran- 
sitional flow states whose temporal complexity changes in 
a rather non-monotonic fashion before a turbulent flow 
is eventually established. 

Kolmogorov flow profile describes a laminar solution of 
the vorticity equation ([2| with the symmetry of the driv- 
ing: continuous translational symmetry in the x direction 
and discrete translational symmetry in the y direction. 
The problem also possesses two additional discrete sym- 
metries (rotation by 180 degrees about a vertical axis and 
a flip about x (or y) axis combined with the change in the 
sign of vorticity), but these will not play an important 
role in the subsequent discussion. The laminar flow (re- 
ferred to simply as L below) is described by the following 
analytical solution for the vorticity 

= — —TTT sin (3) 

and features straight alternating shear bands which re- 
flect the geometric arrangement of the magnets (see Fig. 
[TJa)). For our choice of parameters, linear stability anal- 
ysis predicts this flow profile to be stable for A < 0.1145 



s~^. This is confirmed by the results of our numerical 
simulations summarized in Fig. |2j which shows all stable 
and unstable solutions that have been computed using a 
Jacobian-free Newton-Krylov solver for ^4 < 1 s~^. 

At A ^ 0.1145 s~^ the laminar flow L loses stability 
through a supercritical pitchfork bifurcation and is re- 
placed with its steady, spatially modulated version. As 
A is increased, the distortion of the shear bands increases 
and they are gradually replaced with a periodic array 
of counter-rotating vortices. This spatially modulated 
shear flow (denoted M and shown in Fig. [ijb)) eventu- 
ally undergoes a supercritical Hopf bifurcation and loses 
stability dit A ^ 0.3750 s-^ 

At this point the first stable, time-periodic solution 
(denoted Pi) appears. Four snapshots of this state at 
different phases of the oscillation are shown in Fig. [3] 
For time-periodic flows it is convenient to represent the 




FIG. 3. Time-periodic flow Pi at A = 0.428 s"^ and (a) t=0, 
(b) t=T/4, (c) t=T/2, (d) t=3T/4 with T = 365.83 s. The 
same color bar as in Fig. p[a) is used here. 
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FIG. 4. The perturbation amplitude e and frequency cji = 
27r/T of the time-periodic flows Pi (a), P2 (b), and P3 (c). 
Only the ranges of A are shown where these flows exist and 
are stable. 



stream function as a perturbation about a steady state 

^{x, y, t) = ^o{x, y) + e^i(x, y, t), (4) 

where the perturbation ^1 has zero time average and 
(||^i||2)t = ||^o||2 (II II2 denotes the 2-norm and ( )t 
denotes the time average). The strength e of the time- 
dependent perturbation as a function of A is shown in 
Fig. |4] along with its frequency uoi . 

As expected for a state created via a Hopf bifurcation, 
the amplitude of oscillation for Pi grows as a square root 
of the distance to the bifurcation point (see Fig. [4]^ a)). 
The frequency of oscillations uoi = 27r/T decreases (and 
the period T increases) monotonically with A until the 
oscillatory state is destroyed as a result of an infinite- 
period bifurcation at A ^ 0.4635 s~^. 

At this point two steady solutions are created, a stable 
node N and a saddle S (shown in Figs, [sja) and (b), 
respectively). The corresponding flows are quite similar 
(a disordered array of four clockwise and four counter- 
clockwise vortices) and possess a relatively low symme- 
try: just like Pi, they are symmetric with respect to a 
shift (x, y) ^ {x^ ^3^/2, y + Ly/2). 

The numerical solution of Q follows the stable branch 

as A increases further until the corresponding steady 
flow again develops an oscillatory instability (also a su- 
percritical Hopf) dX A ^ 0.8125 s~^, giving rise to an- 
other time-periodic flow P2, shown in Fig. [6j The ampli- 
tude and frequency of this flow are shown in Fig. [4]^b). 




FIG. 5. Stable steady flow N (a) and unstable steady flow S 
(b) at A = 0.750 s-^ 



This state is stable in a fairly narrow range of A and, at 
A ^ 0.8180 s~^, P2 undergoes a secondary supercritical 
Hopf bifurcation giving rise to a quasi-periodic flow (de- 
noted QP) which, after another Hopf bifurcation, tran- 
sitions to aperiodic flow around A ^ 0.865 s~^. 

At A ^ 0.8740 s~^, a third stable, time-periodic state 
P3, shown in Fig. [7| is created via a subcritical Hopf bi- 
furcation. The corresponding flow does not respect any 
of the symmetries of the system and is only stable for a 
very narrow range of A before it undergoes a subcriti- 
cal pitchfork bifurcation at A ^ 0.8768 s~^. Its ampli- 
tude and frequency are effectively constant throughout 
its range of stability as Fig. ^c) illustrates. 

Increasing A further, we find another narrow aperiodic 
window before the flow returns to quasi-periodic behavior 
at about A ^ 0.885 s~^. Finally, the flow once again 
becomes aperiodic at A 0.980 s~^. The temporally 
aperiodic (or chaotic) flows we find are weakly turbulent. 

We conclude this section by a discussion of the 




FIG. 6. Time-periodic flow P2 at A = 0.817 s"^ and (a) t=0, 
(b) t=T/4, (c) t=T/2, (d) t=3T/4 with T = 131.76 s. The 
same color bar as in Fig. p[a) is used here. 
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FIG. 7. Time-periodic flow P3 at A = 0.875 s"^ and (a) t=0, 
(b) t=T/4, (c) t=T/2, (d) t=3T/4 with T = 94.39 s. The 
same color bar as in Fig. Isjb) is used here. 



Reynolds number 



Re 
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characterizing the solutions described above. As Fig. |8] 
shows, Re varies linearly with A in different flow regimes. 
The slope is roughly the same for almost all flows, except 
the laminar flow L, for which it is much steeper. Indeed, a 
quick inspection of the vorticity fields shows that, beyond 
L, the flow is dominated by structures oriented at an 
angle ^ ^ 45 degrees to the x direction, so that the slope 



can be estimated as Re /A ^ {k/ s'mO)~ 



-2 ^ 47.5 s\ 



For the laminar flow we find instead Re /A ~ k~^i'~'^ ~ 
190 s^. Both estimates are in reasonable agreement with 
the numerical data presented in Fig. [Sj 



III. MIXING PROPERTIES 
A. Numerical Results 

In order to quantify the transport properties of the 
flow, it would be convenient to use two different metrics: 
(i) the relative size (in this case area) of the mixed region 
and (ii) the rate of mixing. Both metrics are most easily 
evaluated by following the evolution of an initially well- 
localized array of passive tracers. Before continuing with 
the detailed discussion of mixing dynamics, we should 
point out that, while the laminar flow L is expected to 
be the worst mixer and the aperiodic (turbulent) flow 
to be the best, the complicated sequence of transitional 
states observed as A is increased implies that we should 
not expect a monotonic increase for either metric. While 
one would expect both metrics to mirror the spatial and 
temporal complexity of the flow, we find that this corre- 
lation is far from perfect. 

As discussed in the introduction, the dynamics of pas- 
sive tracers ([T]) is formally Hamiltonian, with ^ being the 
Hamiltonian. Time-independent, one-degree-of-freedom 
Hamiltonian systems are always integrable and thus ex- 
hibit regular motion. The tracers follow closed stream 
lines on which ^ is exactly conserved, hence, the initial 
tracer distribution eventually stretches along the stream 




A (s-^) 



FIG. 8. Relationship between the Reynolds number and the 
forcing strength A. Again solid lines denote stable states and 
dashed lines denote unstable ones. 



line passing through its center, but never broadens. How- 
ever, the introduction of time-dependence is expected to 
split the flow domain into regions of chaotic and regu- 
lar dynamics. The relation between mixing and chaotic 
stream lines establishes a direct analogy between trans- 
port in one-degree-of-freedom Hamiltonian systems and 
mixing in 2D area-preserving flows. 

In order to quantify the mixing process, for each value 
of A, a set of passive tracers was initially placed in a 
square region with the side of 0.1 mm (which corresponds 
to initial area fraction /(O) = 2 x 10~^). Since the great- 
est degree of stretching usually occurs along homo- or 
heteroclinic trajectories, the initial sets were centered on 
top of one of the saddles of the instantaneous flow fleld. 

Each tracer was then advected by numerically inte- 
grating ([1]) using a fourth-order, area-preserving, sym- 
plectic integrator based on the 2-stage Gauss-Legendre 
scheme [32 . Velocities for each tracer were computed at 
each time step using a cubic interpolation scheme on the 
64 X 128 grid in real space. 

The dispersion of tracers was then used to compute 
the mixing metrics. The mixed area fraction f(t) was 
computed by partitioning the flow domain into a set of 
small boxes and computing the ratio between the number 
of boxes m containing at least one tracer to the total 
number of boxes k. When the tracers uniformly cover 
the domain, the area fraction should be unity. However, 
if there are k boxes with n randomly distributed tracers, 
the fraction of boxes containing at least one tracer would 
on average he Pn,k = 1 — exp(— n//c). Thus, the measured 
area fraction for each value of A was normalized by Pn^k 



fit) 



m{n^ t) 

kpn,k 



(6) 



SO that a uniformly distributed set of tracers would give 
an area fraction of one. 

Fig. [9] shows the area fraction occupied by the tracers 
after a rather long time interval of 5 x 10^ s. In compar- 
ison, the period of Pi, P2 and P3 is of order 100 s, while 
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FIG. 9. The fraction / of the mixed area relative to the total 
area of the domain at t = 5 x 10^ s. 



the characteristic time scale of the flow around vortices is 
below 10 s. We flnd that the area fraction remains near 
zero for all of the time-independent flows (L, M, and N)^ 
as it should be for integrable flows. 

For time-dependent flows ([T]) formally becomes a three- 
dimensional dynamical system (augmented by an equa- 
tion i = 1) which, in general, possesses chaotic solutions 
(stream lines). Chaotic advection, in principle, should 
dramatically enhance mixing. However, as Fig. [9] shows, 
the mixed area fraction for Pi and P2 is only slightly 
higher than that for the time-independent flows. The 
time-periodic flow P3, on the other hand, produces nearly 
perfect mixing, with mixed area fraction comparable to 
that of aperiodic flows. 

Examining the temporal evolution of the area fraction 
covered by the tracers shown in Fig. [lOj one can discern 
two distinct stages for the time periodic flows. Initially 
there is a very fast increase. For Pi and P2 it corresponds 
to rapid stretching of the set of tracers along the homo- 
clinic trajectories forming a thin closed band (see Figs. 
11 'a) and (c)). This is followed by a much slower growth 
associated with the broadening of this band. However, 
even after a very long time, the band of tracers remains 
quite thin and aligned along the stream lines of the in- 
stantaneous flow (see Figs. pTj^b) and (d)). 

For P3, on the other hand, the set of tracers undergoes 
a rapid initial phase of both stretching and folding and 
quickly (within several periods of the flow) covers almost 



the entire domain (see Fig. 
shows that, for P2 and P3, i. 



12[a)). A closer look further 
le tracer distribution reaches 



an asymptotic state already around 10 s, while for Pi 
the area fraction is still growing at t = 5 x 10^ s. Fi- 
nally, although the asymptotic distribution of the tracers 
for P3 is essentially uniform, the tracers never penetrate 
four small regular islands centered around vortices with 
positive vorticity, as Fig , p^b ) illustrates. We will re- 
turn to this fact in Sect. IIII CI 

Fig. [13] shows the tracer distribution for two values 
of A above the onset of the secondary Hopf bifurcation 
which destroys P2 and makes the flow quasi-periodic. We 
flnd the evolution of the tracers to follow the same sce- 
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FIG. 10. Temporal dependence of the area fraction for the 
three time-periodic flows: (a) Pi at A = 0.428 s~^, (b) P2 at 
A = 0.817 s-^ and (c) P3 at A = 0.875 8"^ 



nario as in the case of the time-periodic flow P3: after 
a short initial stage of stretching and folding, the set of 
tracers fllls a signiflcant fraction of the full domain. This 
stage is followed by a much slower homogenization pro- 
cess in which the distribution becomes spatially uniform. 
However, just like in the case of P3, the tracers never 
penetrate four regular islands centered around vortices, 
now with negative vorticity. 

The fundamental difference between (quasi) periodic 
and aperiodic flows makes itself apparent if we compare 
mixing by the periodic flow P3 with that by aperiodic 
flows just outside of the window of stability for P3, at 
A = 0.872 s-2 and A = 0.878 s-^ Although the forcing 
is almost identical in these three cases and the short- 
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term dynamics of the three flows are similar. Fig. 
shows that the aperiodic flows achieve perfect mixing in 
the long term, covering the entire domain, including the 
four regular islands of P3. 



Fig. 15 summarizes the observed mixing rate as a func- 
tion of the control parameter A. As Fig. [To] amply illus- 
trates, the mixing process is characterized by a range of 
time scales. The fastest time scale describes stretching of 
the initial tracer distribution along the stream line pass- 
ing through its center. The corresponding rate is deflned 
as Tmax = max^ \df/dt\ and is proportional to the average 
shear rate corresponding to that stream line. 

The slowest time scale describes broadening of the 
distribution due to transport of tracers through semi- 



penetrable transport barriers discussed in Sect. HI C To 
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FIG. 11. Mixing by the time periodic flows. The distribution of 6 x 10^ tracers and the stream hnes of the instantaneous flow 
for Pi at t = 598 s (a) and t = 5 x 10^ s (b). The same for Ps at t = 49 s (c) and t = 5 x 10^ s (d). 
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FIG. 12. Mixing by the time-periodic flow P3. The distribution of 6 x 10^ tracers and the stream hnes of the instantaneous 
flow at t = 317 s (a) and t = 3500 s (b). 



characterize this broadening, we computed the time tgo 
it takes for the area fraction to reach 90% of its asymp- 
totic value, f (too) / f {hoo) = 0.9, where we assumed the 
asymptotic distribution is achieved at tioo = 5 x 10^ s. 
The minimal mixing rate was then defined as rmin = 
l/tgo- In both cases we averaged f{t) over a small win- 
dow to filter out small oscillations associated with the 
passage of tracers near saddles. 

The fast time scale rmax is found to increase almost 
monotonically with A, reflecting the corresponding in- 
crease in the shear of the underlying flow. The slow time 
scale requires more care to interpret. In particular, for 
Pi we find rmin to drop by almost an order of magnitude 
as A increases. This decline is associated with the tracer 
distribution shown in Fig. [lllb) slowly broadening in 
time as illustrated by Fig. |lQ[a). This broadening is 
due to a slow "leak" of tracers across a semi-penetrable 
transport barrier, creating a "halo" of tracers surround- 



ing the main band. Another drop observed around the 
secondary Hopf bifurcation at A ^ 0.818 s~^ is associ- 
ated with a similar process for the quasi-periodic flow 
QP. As A increases past this critical value of A, the 
transport barrier which exists for P2 gets eroded, leading 
to a quick increase in rmin- 

While many of our numerical results are quite logical, 
several findings raise questions. For instance, the flows Pi 
and P3 appear to be qualitatively quite similar. Both are 
stable, time-periodic and, with the choice of A = 0.428 
s~^ for Pi, both have a time-dependent component of 
the same magnitude e ^ 0.238. Yet, despite these sim- 
ilarities, their mixing properties are radically different. 
Pi is a very poor mixer, as Fig. pTj^b) illustrates. It is 
characterized by both a very low mixing rate and a very 
low mixed area fraction. In fact. Pi's mixing properties 
are comparable to those of time-independent flows. P3, 
on the other hand, is an extremely good mixer, almost 
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as good as the aperiodic flows. The mixing rate for this 
flow is high and its mixed area fraction is close to unity. 

Another question concerns the islands surrounding 
positive or negative vortices that remain impenetrable 
for extremely long times for both the time-periodic flow 
P3 (Fig. and the quasi-periodic flow QP succeeding 
P2 (Fig. . In both cases there appear to be transport 
barriers surrounding vortices characterized by vorticity 
of one sign but not the other. This was also found to 
occur in the model flow of Danilov et al. [25 as well as 
in real oceanic flows l33l. 



B. Lagrangian Coherent Structures 

Finite-time Lyapunov exponent (FTLE) flelds associ- 
ated with the time-dependent flows provide some intu- 
ition regarding their drastically different mixing proper- 
ties. The forward FTLE is a scalar quantity 



cr(xo,r) = - In 
r 



Dx(r) 



(7) 



which characterizes the amount of stretching along a tra- 
jectory x(t) passing through the point xq at t = over a 
flnite time interval r. In particular, the ridges of the for- 
ward FTLE fleld deflne Lagrangian Coherent Structures 
(LCS) [34] which, for time-periodic flows, correspond to 
segments of unstable manifolds of saddle orbits with tem- 



poral period equal to that of the flow. As Fig. [16] illus- 
trates, for Pi and P2 the LCS show very little folding, 
effectively forming closed, compact curves. For P3, on 



the other hand, the LCS display a lot of folding, which 
is a necessary ingredient for efficient mixing and cover a 
substantial fraction of the total area. Indeed, we flnd the 
LCS of Pi and P2 are qualitatively similar to those of 
steady flows from which they are born (i.e., M and A/"), 
while the LCS of P3 are qualitatively similar to those 
of aperiodic flows, which is consistent with the observed 
similarities in their mixing properties. 

LCS play an important role in organizing transport. 
For instance, placing the initial set of tracers on top 
of the saddle orbit we should expect that set to be 
quickly stretched along the LCS forming effectively one- 
dimensional structures for Pi and P2, while for P3 the 
structure becomes effectively two-dimensional. Further- 
more, the LCS form transport barriers which cannot be 
crossed by the tracers. For Pi and P2 (as well for steady 
flows), these transport barriers are closed, effectively par- 
titioning the domain and preventing mixing between re- 
gions separated by the LCS. For P3 (as well as for ape- 
riodic flows), the transport barriers are open, enabling 
transport and mixing across the whole domain. 

The LCS-based description of transport is consistent 
with our long-term numerical advection calculations and 
has the advantage that it requires time-integration over a 
considerably shorter time-interval (fraction of the tempo- 
ral period T of the flow, compared with hundreds to thou- 
sands of periods for numerical advection calculations). 
However, neither approach explains why the mixing prop- 
erties of the time-periodic flows are so dramatically dif- 
ferent. A more insightful approach is discussed next. 



9 





C. Resonance Phenomena 

As we mentioned previously, area-preserving time- 
periodic flows Pi, P2, and P3 can be treated formally as a 
perturbed Hamiltonian system ([T]) , with the stream func- 
tion Q serving the role of the Hamiltonian. In particu- 
lar, ^0 plays the role of the unperturbed Hamiltonian and 
e^i - the time-periodic perturbation. Transport in near- 
integrable time-periodic Hamiltonian systems and area- 
preserving flows has been studied extensively. It is well 
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FIG. 15. The rates of mixing for time-dependent flows as 
a function of A. The solid and dashed curves correspond, 
respectively, the fastest time scale rmax and the slowest time 
scale rmin- 



understood that, for weak perturbations, chaotic trajec- 
tories emerge in the neighborhood of the homo- or hetero- 
clinic manifolds of saddle orbits of the integrable unper- 
turbed, or base, flow. These manifolds self-intersect as a 
result of the imposed perturbation, forming a homoclinic 
tangle with the lobe dynamics [19^ which provides an in- 
sightful, albeit computationally challenging, description 
of mixing in the separatrix chaotic layer (SCL). 

Not only is the computation of the width of the SCL 
is a computationally intractable problem for any real- 
istic flow, the width of the SCL significantly underes- 
timates the size of the actual chaotic domain for finite 
values of e. According to the KAM theory [35-37 , in 
the presence of perturbation, resonant tori of the unper- 
turbed flow (tori whose frequency cJo(^o) is in rational 
ratio with the frequency of the perturbation uji = ^ir/T) 
break up, forming chains of elliptic and hyperbolic time- 
periodic orbits (or stream lines) with their own sets of 
self-intersecting stable and unstable manifolds generating 
resonant chaotic layers (RCL). These RCLs can overlap 
with the SCL making the chaotic domain much broader. 

The dynamics away from the separatrix can be de- 
scribed by computing the change in the value of ^0 over 
an interval of time (0,t/). By analogy with the deriva- 
tion of Melnikov's function, we can use Q and ([T]) to 
show that 

^o(^/)-^o(0)=6 [ \o{^{t)W{^{t),t)dt, (8) 

where = {dy^i^ —dx^i) and the superscript ± denotes 
the component normal to the stream line of the unper- 
turbed flow. The velocity field describing a time-periodic 
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FIG. 16. Forward finite-time Lyapunov exponent field, (a) 
Pi at A = 0.428 s"^ witii r = 32 s, (b) P2 at A = 0.817 s"^ 
witii r = 19 s, and (c) P3 at A = 0.875 s"^ witii r = 22 s. 



perturbation can be written in the form of a Fourier series 
vi{^,t) = J29i^,^k)e-'^'*, (9) 

k 

where we have defined ujk = kuj\. Furthermore, the prod- 
uct 'Uo(x(t))^(x,a;/e) is also a time-periodic function with 
period Tq = 2ti /uoq and can be written as a Fourier series 



imuJot 



(10) 



Substituting and (10) into (|8| we find that, over a 
time interval tj ^ max(T, To ) , the rate of change 



*o(i/) -*o(0) 



vanishes unless ojq/oji = k/m for some integer m and k. 
In what follows, we take m and k to be positive. 

The Fourier coefficient Gk,m also controls the width 
Wl ^ of the RCL which replaces the stream line of the 



0.5 r 



3 0.25 





100 

t (s) 



200 



FIG. 17. Velocity magnitude for a typical stream line near a 
separatrix of the base flow of Pi . 



unperturbed flow with frequency uJo = ook/m. For stream 
lines close to the separatrix '^o{x,y) = tpi^ 'Uo(x(t)) is 
small everywhere except for short time intervals corre- 
sponding to fast motion away from the saddles (see Fig. 
17), so we can estimate 



\G, 



k,m I 



1 

- / vo(x(t))^(x(t),cj/e)e 
^0 Jo 



-imujQt 



dt 



si ujk\g(yi],uJk)\ 



27r m 
where si is the length of the separatrix. 



vii^*,t)e 



(12) 



(13) 



is the (discrete) Fourier spectrum of 'U^(x*,t), and x^* is 
the point on the separatrix which lies midway between 
the saddles (for which i'o(x(t)) is near its maximum). 
If there is more than one saddle on the separatrix, the 
estimate ([l2| should be generalized to include respec- 
tive contributions from all segments, which can either 
enhance or suppress each other. 

The period of the unperturbed motion along the 
stream line lying near the same separatrix is given by 



To(*o) = -^A-/ln 



6 



(14) 



where A^^^ are the positive eigenvalues of all the saddles on 
the separatrix and is a constant. Hence, the distance 
(in terms of ^o) from the separatrix to the nearest k'.m 
resonant torus is exponentially small for low k: 



ii exp 



ii exp 



vfvxx\ 
kuJi J ' 



(15) 



where xT^ = ^i^l l"^^- ^ similar fashion we can 
compute the distance between various resonant tori. In 
particular, the distance between the tori with frequency 
ratios k:l and (A: + 1):1 is given by 



Si{uJk) 



dujo 



exp 



_Xi_ 
^k 



(16) 



According to (12), for moderate /c, \Gk^ra\ takes the 
largest values for m = 1, hence the width of the dominant 
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FIG. 18. Frequency of the motion along the stream hnes of the base flow for (a) Pi at A = 0.428 s~ , (b) P2 at A = 0.817 
s~^, (c) P3 at A = 0.875 s~^ and (d) QP at A = 0.820 s~^. The black portion corresponds to the chaotic domain around 
the separatrix ^0 = seeded with the tracers, while the gray portion corresponds to the regular (as well as chaotic) domains 
without tracers. Horizontal lines show the frequencies of the overlapping dominant RCLs. 



{k:l) RCLs can be estimated from (11): 



''4\Gk,i\ 



esi uJk 
2 uji 



\g{^huj,)\. (17) 



Comparing the widths of the RCLs with the distances 
Si{uJk) between the neighboring resonant tori we can de- 
termine which RCLs overlap and which do not for a par- 
ticular strength of the perturbation. 

For moderate e, we can expect several RCLs with low k 
to overlap with each other and with the SCL, since Si{uJk) 
is exponentially small, while Wi{uJk) scales as a power 
of UJk near a separatrix. More specifically, the region 
where Wi{uJk) > Si{uJk) is expected to be well mixed, 
while in the region where Wi{uJk) < Si{uJk) mixing is 
expected be limited to narrow RCLs of width Wi{uJk) (as 
well as some even narrower RCLs corresponding to m > 
1). The boundaries of the main chaotic region should be 
determined by the regular tori with frequency 



duJn 



The value k± on each side is different and corresponds to 
the outermost overlapping RCL, i.e., it is largest integer 
k± such that Si{uJk) < Wi{uJk) for all k < k±. 



Consider, for example, the flow P2. Fig. 19 ^b) shows 
the width of resonant chaotic layers and the spacing be- 
tween the resonant tori on both sides of the saddle seeded 
with tracers (which corresponds to ^0 = 0). For low 
values of k we do indeed find Wi{uJk) > Si{uJk), so the 
separatrix chaotic layer and the RCLs of a few nearby 
k:l resonant tori overlap, forming a single chaotic do- 
main. Wi{uJk) decreases quickly (in fact, exponentially 
fast) with cj/e, while Si{uJk) increases on both sides of the 
separatrix, so that the number of overlapping RCLs is 
rather low (/c_ = 3 for ^0 < and /c+ = 4 for ^ > 0), 
making the width of the chaotic domain (the black region 
in Fig. [Tsj^b)) extremely small, in good agreement with 
the numerical result shown in Fig. [lljd). 



The situation is similar for Pi. Fig. |l9[a) shows that 
the number of overlapping RCLs is somewhat larger for 
Pi {k- = k-^ = 11), but still small enough for the chaotic 
domain formed by the overlapping resonant and separa- 
trix chaotic layers to remain quite thin (see Fig. [isj^a)). 
This is again in agreement with the numerical result 

n;b)- 



(18) shown in Fig. 



11 



A "halo" of tracers outside of the 
well mixed region suggests that the outermost RCL just 
touches its neighbor on the side of the separatrix, with 
the two separated by a semi-penetrable transport bar- 
rier formed by either a narrow high-order RCL with very 
small Gk,m or by a cantorus [38j, which allows a very 
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FIG. 19. The widths Wi{uj) of resonant chaotic layers (soUd Une) and the spacings Si{uj) between the corresponding resonant 
tori (dashed and dotted hnes) on different sides of the separatrix = 0. (a) Pi at A = 0.428 s~^ with uji = 0.0172 s~^, (b) 
P2 at A = 0.817 s"^ with uji = 0.0477 s"\ (c) P3 at A = 0.875 s"^ with uji = 0.0666 s"\ (d) QP at A = 0.820 s"^ with 
cji ^ 0.051 s"^ 



slow "leak" of tracers into the outermost RCL. 

While Pi and P2 are almost monochromatic, the 
Fourier spectrum ^(x*,^;^) of P3 is exceptionally broad, 
with a large number of harmonics cuk that have ampli- 
tudes comparable to that of the base frequency uji. As 
a result, we find that RCLs remain fairly broad even for 
large values of k (see Fig. [l9|c)). Consequently, the 
chaotic domain for P3 is comprised of a large number of 
overlapping RCLs which allow transport across most of 
the ^0 range (that is across most of the physical space). 
Since in this case the boundaries of the chaotic domain 
are very far from the separatrix of the saddle seeded with 
tracers, we computed the spacing Si{uJk) between the res- 
onant tori corresponding to the two outermost branches 
of the co'o(^o) curve which extend to the extremal val- 
ues of ^0 correspondin^to clockwise and counterclock- 
wise vortices. As Fig. 
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'c) illustrates, Wi{uJk) > Si{uJk) 
for all k for both the leftmost and the rightmost branch. 
Hence, we should expect all of the RCLs to overlap, allow- 
ing global transport. However, the last RCL surround- 
ing the clockwise vortex is not wide enough to cover the 
whole range of leaving a small regular island around 
each of the four vortices centered at '^o{x,y) = 1.058, 
which is also in agreement with the numerical result 
shown in Fig. 12 'b). 

The mixing properties of the quasi-periodic flow QP 
can also be understood by analyzing the widths of RCLs. 



The time-average flow for QP is essentially the same as 
that for P2, hence we can use the same frequency curve 
a;o(^o)- The spectrum ^(x*,cj) of the quasi-periodic per- 
turbation is discrete, just like the spectra of the periodic 
flows, with frequencies of the peaks that can be labeled 
cj/e, with integer k. Although the spacing cj^+i — be- 
tween the peaks is not exactly constant, it varies little 
about the average uji ~ 0.051 s~^, as Fig. 19 'd) illus- 
trates. Comparison of Si{uJk) and Wi{uJk) shows that 
RCLs with k < 9 overlap. In addition, the tori 13:1, 14:1, 
15:1, and 16:1 also overlap. Although the width of the 
tori 10:1, 11:1, and 12:1 is smaller than the spacing be- 
tween them, the tori 9:1 and 13:1 are more than twice as 
wide as the biggest inter-tori spacing in this "gap" , which 
means that the 9:1 and 13:1 RCLs overlap directly. This 
indicates that we should have global transport in the re- 
gion of ^0 where a;o(^o) < ^16 + Wi{uJiQ)\duJo/d'^o\/2 ^ 
0.85. According to Fig. [iSjd), this corresponds to almost 
the entire physical domain, with the exception of regular 
islands around all the vortices. 

Although this prediction is not in perfect agreement 
with the numerical result, the description in terms of 
interacting resonances captures most of the features of 
the asymptotic tracer distribution shown in Fig. ^h). 
Both global transport and the regular islands around the 
clockwise vortices are predicted correctly. The size of the 
regular islands is predicted to be much larger than for 
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Pa? which is also consistent with numerics. The resonant 
description also predicts the formation of regular islands 
around the counterclockwise vortices, which are filled in 
in Fig. 13, which shows a limitation of our analytical 



description. 

However, the discrepancies at high frequencies are ex- 
pected, given the fact that the expression (12) for the 
Fourier coefficients Gfc,m (and hence the widths Wi{uJk)) 
was obtained in the limit of low frequencies cj^. For 
higher frequencies corresponding to the neighborhood of 
the vortices, ([T2| becomes inaccurate and Gk,m has to be 
computed in a different way. 



IV. SUMMARY AND CONCLUSION 

To summarize, we have described transition from the 
laminar Kolmogorov flow to turbulence in a doubly- 
periodic domain of relatively small size. The sequence 
of bifurcations preceding turbulence is quite rich, with 
several different steady and time-periodic flows succeed- 
ing one another. This bifurcation sequence is quite sen- 
sitive to the choice of parameters (a, /3, the system 
size Lx X Ly) although the actual transition to turbu- 
lence follows one of two standard routes. In one case 
we find that turbulence emerges through a sequence of 
several Hopf bifurcations, commonly referred to as the 
Ruelle-Takens-Newhouse scenario [W, "W. In the other, 
a sub-critical bifurcation leads to intermittency, as in the 
Pomeau-Mannevile scenario [41]. 

The details of the bifurcation sequence, however, are 
quite important in describing the evolution of the trans- 
port properties of the flow. As a general trend, we find 
the mixing efficiency (defined either in terms of the mixed 



area fraction or in terms of the mixing rate) to improve 
as the forcing is increased, with steady flows being the 
worst mixers and turbulent flows - the best. However, 
the complexity of the flow does not increase monotoni- 
cally and neither does mixing efficiency. Neither is mix- 
ing efficiency directly related to the complexity of the 
flow, as the comparison of three different time-periodic 
flows showed. Furthermore, time-periodic flows such as 
P3 can rival the mixing efficiency of turbulent flows. 

The most unexpected result was that the mixed area 
fraction of a class of time-periodic and quasi-periodic flow 
can be described - rather accurately - by a perturbative 
approach. This is despite the fact that none of the flows 
considered can actually be considered weakly perturbed. 
The description is based on the idea of multiple over- 
lapping resonances, which define well-mixed regions of 
the flow. In particular, our results conffim the idea of 
Soskin and Mannella [42 that resonances play an impor- 
tant role in defining the width of (and dynamics inside) 
the separatrix chaotic layer. Although the flow domain 
is densely covered by an infinite number of resonant tori, 
only the dominant resonances (e.g., ones that correspond 
to the harmonics of the frequency of the perturbation) 
play an important role. As a general rule of thumb, we 
find the flows with the broadest Fourier spectrum to pos- 
sess the best mixing properties, while (nearly) monochro- 
matic flows have mixing properties comparable to those 
of steady flows. 
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